Fast computation of multi-scale combustion systems 
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In the present work, we illustrate the process of constructing a simplified model for complex 
multi-scale combustion systems. To this end, reduced models of homogeneous ideal gas mixtures 
of methane and air are first obtained by the novel Relaxation Redistribution Method (RRM) and 
thereafter used for the extraction of all the missing variables in a reactive flow simulation with a 
global reaction model. 
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INTRODUCTION AND MOTIVATION 

Solution of the full set of equations as required in nu- 
merical simulations of reactive flows with detailed chem- 
ical kinetics represents a quite challenging task even for 
today super-computers. The one reason is the large num- 
ber of kinetic equations needed for tracking each chemical 
species. On the other side, detailed combustion mech- 
anisms are typical multi-scale problems where different 
chemical processes, characterized by disparate timescales 
ranging over several orders of magnitude (from seconds 
down to nanoseconds) are present. As a result, model- 
ing detailed combustion fields comes with a tremendous 
cost where intensive long simulations are needed to re- 
solve the fastest processes, though one is often interested 
in the slow dynamics. Thus simplification methodologies 
become, other than highly desirable, mandatory in com- 
bustion problems where detailed mechanisms for heavy 
hydrocarbons (with hundreds chemical species) are used 
in 2- and 3D simulations. 

Notice that, there exist often chemical processes that 
are much faster than the fluid dynamic phenomena, so if 
we are only interested in computing the system behavior 
on the time-scale of the fluid mechanics, some chemical 
processes will be already equilibrated and thus slaved to 
the remaining dynamics. 

In fact, modern simplification techniques are based on 
a systematic decoupling of the fast equilibrating chemical 
processes from the rest of the dynamics, and are typically 
implemented by seeking a low dimensional manifold of 
slow motions in the solution space of the detailed system. 

Much effort has been devoted to setting up such au- 
tomated model reduction procedures. The method of 
mvar iant grids (MIG) 0, Q, the computational singu- 
lar perturbation (CSP) method Q, the intrinsic low di- 
mensional manifold (ILDM ) Q , the invariant constrained 
equilibrium edge preimage curve method (ICE-PIC) fioj ]. 
and the method of minimal entropy production trajecto- 
ries (MEPT) are a few popular techniques. 

In this work, we introduce an approximated procedure 
for the fast computation of detailed combustion fields. 



To this end we adopt the novel Relaxation Redistribu- 
tion Method (RRM) for the construction of a reduced 
model of the mechanism GriMech 3.0 describing ideal 
mixtures of air and methane (53 chemical species, 4 ele- 
ments and 325 reactions) in a closed system under fixed 
pressure and mixture-averaged enthalpy. The latter de- 
scription serves, at a later time, for reconstructing all the 
missing chemical species in a computationally efficient 
reactive flow simulation performed with a single-step re- 
action model. 

The paper is organized in sections as follows. The no- 
tion of slow invariant manifold and the chemical kinetics 
equations are briefly reviewed in section „ and „ respec- 
tively. The relaxation redistribution method is discussed 
in section n A reduced model for air and methane is used 
in a planar counter-flow flame simulation in section,, and 
conclusions are drawn in section „ 



SLOW INVARIANT MANIFOLD (SIM) 

In this section, we briefly discuss the notions of slow 
invariant manifold for a system of autonomous ordinary 
differential equations in a domain 14 in BP , 



V = f(y)- 



(1) 



For more details, the interested reader is delegated to the 
dedicated literature 0, Q . A manifold Q C U is invariant 
with respect to the system ([I]) if inclusion y(to) G O 
implies that y(t) £f! for all future time t > to- 

Equivalcntly, if the tangent space T y to f2 is defined at 
y, invariance requires: f(y) G T y . In order to transform 
the latter condition into an equation, it proves conve- 
nient to introduce projector operators. Let for any sub- 
space T y a projector P onto T y be defined with image 
imP = T y . Then the necessary differential condition can 
be expressed by: 



(l-P)/ = 0, 



(2) 



where the left-hand side of equation [2] is often called 
defect of invariance A. It is worth stressing that, al- 
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FIG. 1: Basic idea behind the Relaxation Redistribution 
Method. 



though the notion of invariance discussed above is rel- 
atively straightforward, slowness instead is much more 
delicate. We just notice that the most part of invariant 
manifolds are not suitable for model reduction (all semi- 
trajectory are, by a definition, ID invariant manifold). 
In this respect, we should also point out that slow man- 
ifolds are not uniquely defined in the literature and, in 
general, different methods delivers different objects. 

Here, we follow the approach of the Method of Invari- 
ant Grid (MIG) where slowness is intended as sta- 
bility, thus SIM are the stable stationary solution of a 
relaxation process (film equation) 



dF(Q 
dt 



(1 - P)f- 



(3) 



We notice that the unknown of both ^ and Q is the 
manifold fi, which can be conveniently represented in a 
parametric form, as mapping F : W — > U of a domain W 
(reduced space or parameter space in the following) into 
the phase-space IA; f2 being the image of this mapping: 
O = F(W). 



REACTION KINETICS EQUATIONS 



and specific enthalpy of species i, respectively. Following 
ChemKin [B[ , the specific enthalpy can be approximated 
by a polynomial fit as follows: 



hi(T) = RT(a u 
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(5) 

where aji are the tabulated Nasa coefficient and R the 
universal gas constant. Molar concentration rates take 
the explicit form: 



3 = 1 



(0) 



where oiij and ftj are the stoichiometric coefficients of 
the j-th elementary reaction X)"=i a ijXi =± Xa=i 
The j-th reaction rate rj is expressed using the popular 
mass action law. 



k f n - k~ ii [x, 



(7) 



with [Xi] denoting the molar concentration of the ith 
species and k^ the j-th reaction rate constant typically 
expressed in the Arrhenius form [13| . 



MODEL REDUCTION TECHNIQUE 

In the following, we use a discrete representation for 
manifolds referred to as grids Q , consisting of a set of in- 
terconnected nodes, where it is assumed that the nearest 
neighbors of an arbitrary node y can be identified. A grid 
is defined by the restriction of mapping F on the discrete 
subset of the parameter space GcW into the phase space 
IA, whereas and invariant grid satisfies the grid version of 
the invariance equation: / \F '(£)) - P f (F '(£)) = 0, V£ e G 
jij. Notice that, thanks to the node connectivity, it is pos- 
sible to compute local tangent space hence the projector 
P (e.g. using approximated differential operators). 



Here, we consider closed reactive systems with fixed 
mixture-averaged enthalpy h, and total pressure p where 
n chemical species and d elements participate in a com- 
plex network of elementary reactions. Species compo- 
sitions are represented in terms of mass fractions Yi = 
mi/mtot, with mi and m to t denoting the mass of species 
i and the total mass, respectively. The mixture en- 
thalpy, at a temperature T, can be expressed as h = 
J27=i h-iiT^i, while t he g overning equations in a closed 
reactor take the form [13j : 



y = f(y) = ( 



LJiWi 
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LO n W n T 



(4) 



where p is the mixture density, while Wi, LOi, hi(T) de- 
note the molecular weight the molar concentration rate 



Relaxation methods 

Here, construction of one-dimensional invariant grids is 
accomplished by the Relaxation Redistribution Method 
(RRM) which has proven an efficient method for solving 
the film equation ((3]) starting from an initial grid Qq. The 
interested reader can find further details in 1]. Referring 
to Fig. [TJ for simplicity, in this work C?o is chosen regular 
in terms of the parameter £. 

Let a numerical scheme (Euler, Runge-Kutta, etc.) be 
chosen for solving the system of kinetic equations (|4|), 
and let all the grid nodes relax towards the slow invari- 
ant manifold (SIM) under the detailed dynamics / dur- 
ing one time step. The fast component of / brings a grid 
node closer to the SIM while, at the same time, the slow 
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FIG. 2: Relaxation redistribution method where an explicit 
first order Euler scheme is adopted during relaxation with an 
adaptive time step St < 1.5 x 10 -8 . 



component causes a contraction towards the steady state 
of (|4]). As a result, the grid becomes dense in a neigh- 
borhood of the steady state and coarse far from it, when 
keeping relaxing. The slow motion can be neutralized 
by a node redistribution after the grid relaxation (thus 
mimicking the term —Pf in ([3])). In other words, as il- 
lustrated in Fig. [TJ the relaxed states are redistributed 
on a regular grid in terms of the parameter £ via linear 
interpolation. 

Notice that, all intermediate grids are, by construction, 
regular in terms of £ and, in the case of an invariant grid, 
the overall effect due to relaxation and redistribution is 
null and the invariance condition satisfied. 

In our study, we consider a detailed combustion mech- 
anism for air and methane (GriMech 3.0), where n = 53 
chemical species and 4 elements participate in 325 ele- 
mentary reactions Here, at a fixed mixture-averaged 
enthalpy h and pressure p, Go represents the mixing line 
between the two states yf resh (stoichiometric fresh mix- 
ture) and y eq (stoichiometric chemical equilibrium state) 
discretized by N = 200 nodes. Iterations have carried 
out until Sy/Sy < e at every grid node, where Sy is the 
overall movement due to the relaxation and redistribu- 
tion while Sy is the movement due to relaxation alone of 
an arbitrary node y with a tolerance s = 0.001. To the 
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FIG. 3: Yco coordinate. Two-dimensional invariant grid via 
the relaxation redistribution method. 



end of constructing a two-dimensional invariant grid pa- 
rameterized with respect to £ = Ycoi and h, the above 



construction is performed over a range of enthalpies —415 
[kJ/kg] < h < -5 [kJ/kg] with a step Ah = -5 [kJ/kg]. 
In Figurc|3l a projection of the above invariant grid in the 
three-dimensional sub-space h — Yqoi ~ Yco is reported. 



REACTIVE FLOW SIMULATION 
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FIG. 4: Planar counter-flow configuration. Symmetry is as- 
sumed with respect to the line y — L. 

Let us consider the planar stagnation point flow, where 
a well premixed stoichiometric mixture of air and fuel, 
initially at room condition (T = 300A', p = lbar), im- 
pinges against a stream of hot products. Due to sym- 
metry, a flat flame can be established in this flow at 
< y < L as schematically depicted in Figure [4] Al- 
though the above is effectively a two dimensional prob- 
lem, under the assumptions of symmetry, boundary layer 
approximation and low Mach number regime, it is possi- 
ble to consider the following one dimensional system of 
governing equations imposing conservation of mass, mo- 
mentum, energy and chemical species, respectively, along 
the symmetry- line (x = 0): 



dp dV_ 

dt dy 



pUe = 0, 



dU Tr2 , r dU 



dy dy 



pfresh e 



(8) 



0, (9) 



_dT dT u u± ^ 

P ~dt + ~dy ~ T^-Th,^ Jh,> 2^ 



1__5 

C p dy K 



dT, 
'dy 1 



LUihi 



= 0, 



(10) 



.dY, 8Yi 

P^T- + V- 



dy {pDi ^ ] 



WiWi = 0, 



l,...,n, 



dt dy 

(11) 

where the ideal gas law, p = pW/RT, can be used for the 
closure. Moreover, p, A, Di and pf resh are the dynamic 
viscosity, thermal conductivity, the diffusion coefficient of 
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species i and the density of the fresh mixture at the inlet, 
respectively. The mean specific heat C p (under constant 
pressure) and the mean molecular weight W take the 
explicit form: 

" - 1 

C p = W = Eti y t/Wi - 

with c P i being the specific heat of species i (mass unit). 
Let u, v and e be the two velocity components of 2D flow 
field along the x and y axes and the flame strain rate, 
respectively. We note that the above set of equations ((8])- 
are conveniently expressed in terms of the quantities 
U = u/uoo and V = pv, with — ex. The latter 
problem is solved imposing fresh mixture condition at 
the inlet (y = 0): 

y = 0: U = l, T = 300[K], Y i = Y i lreah , (12) 

while chemical equilibrium and zero-flux condition can 
be chosen at the outlet 

V = L: = 0,^ = 0,^ = 0, 

V P eq oy oy 

(13) 

where p eq is the density of the fully burned mixture. 

The detailed derivation of the set of equations (j5])- (fTTj) 
along with the boundary conditions (fT2")l and (fT3"|) can be 
found in Q. In this study, spatial derivatives in (j51)- (fTTj) 
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FIG. 5: Premixed flame in a planar counter-flow. Equa- 
tions (f8))- (|lip are solved by finite difference scheme using 
the global reaction model (|14p and (|15p . Minor chemical 
species are thereafter computed accessing the invariant grid 
in Figure O Here, we use A = D = 3 X l(T 5 [m 2 s _1 ], 
A = 0.026[iym -1 ^ -1 ], fixed p = l[bar] and strain rate 
e = 401s" 1 ]. 

have been approximated by finite differences (upwind for 
convective terms while central differences for diffusive 
terms) , and the corresponding ordinary differential equa- 
tions (ODE) system has been solved by a numerical stiff 
solver [llj readily available in Matlab® (odel5s). More- 
over, we first consider n = 4 reactive chemical species 



(CH4, CO2, O2, H2O) with an abundant inert (N2) par- 
ticipating in the one-step global oxidation (s = 1): 

CHi + 20 2 + 7.52N 2 -)■ C0 2 + 2H 2 + 7.527V 2 , (14) 

whose rate r\ , according to [l3| , can be evaluated as fol- 
lows: 

n = -1.3 x io 8 e- 24358 / T [CH 4 ]-°- 3 [0 2 ] 13 . (15) 



Missing fields retrieval 

Upon the solution of the governing equations (|51)- (fTTj) 
in combination with the global reaction mechanism (|14l) , 
five chemical fields, the temperature, the mixture den- 
sity and flow velocity are available along the symmetry 
line in Figure [4] On the one hand, the computational 
cost of the latter simulation is drastically reduced com- 
pared to a reactive simulation with a detailed combustion 
mechanism, where a much stiffcr and larger set of species 
equations (|11[) must be taken into account. However, on 
the other hand, global approaches inevitably come with 
a significant lack of information concerning all intermedi- 
ate and minor chemical species which underlie a complex 
phenomenon such as the one represented by (114[) . 

In order to fill this gap, here we suggest to perform an 
a posteriori retrieval of all the missing variables in the 
above computation (e.g., minor species such as radicals) 
via linear interpolation from the table describing the in- 
variant grid evaluated (once and forever) in section „ by 
means of the detailed combustion mechanism GriMcch 
3.0 [IH {n = 53 and s = 325). To this end, the two vari- 
ables Yco 2 1 h can be extracted from the above simulation 
and used to access any of the coordinates of the invariant 
grid (see also Figure [3]) . Some of the interpolated vari- 
ables, in a well developed flame along the channel, are 
reported in Figure [SJ 



CONCLUSION 

In this work, we first demonstrate that a recently 
introduced model reduction technique (Relaxation Re- 
distribution Method - RRM) is suitable for handling a 
complex multi-scale combustion mechanism for hydro- 
carbons. Moreover, on the basis of the latter simplified 
model, we introduce and test an embarrassingly simple 
method for the computation of minor chemical species 
upon a reactive flow simulation efficiently performed with 
a global (not necessarily one-step) reaction. The present 
study is the first step towards fast computation of de- 
tailed combustion fields for heavy hydrocarbons in 2D 
and 3D problems. 
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